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We present a numerical self consistent variational approach based on the Jordan- Wigner transfor- 
mation for two dimensional spin systems. We apply it to the study of the well known quantum 
(S — 1/2) antiferromagnetic XXZ system as a function of the easy-axis anisotropy A on a pe- 
riodic square lattice. For the SU(2) case the method converges to a Neel ordered ground state 
irrespectively of the input density profile used and in accordance with other studies. This shows the 
potential utility of the proposed method to investigate more complicated situations like frustrated 
or disordered systems. 

I. INTRODUCTION 

Quantum spin systems in two dimensional (2D) lattices have been the subject of intense research, mainly motivated 
by their possible relevance in the study of high temperature superconductors^. On the other hand, high magnetic 
field experiments on materials with a 2D structure which can be described by the Heiscnbcrg antiferromagnetic model 
in frustrated lattices have revealed novel phases as plateaux and jumps in the magnetization curves^. In spite of the 
huge efforts made, a general understanding of the phase diagram of such magnets is elusive and it is then worth trying 
to develop new techniques to study these systems systematically. Among the many different techniques that have 
been used to study such systems, the generalization of the celebrated Jordan- Wigner (JW) transformation^ to two 
spatial dimensions^ has some appealing features. It allows one to write the spin Hamiltonian completely in terms of 
spinlcss fermions in such a way that the S = 1/2 single particle constraint is automatically satisfied due to the Pauli 
principle, while the magnetic field enters as the chemical potential for the JW fermions. The price one has to pay is 
the appearance of complicated non-local interactions between fermions. This method has been applied in£ (see also£) 
to study the XXZ Heiscnbcrg antifcrromagnct. These studies have been reviewed int 

More recently this technique was used to obtain a theoretical magnetization curve for the Shashtry-Sutherland 
model, reproducing at the mean field level some of the experimentally observed features for the material SrCu2 (1303)2 
which is assumed to be described by such model^. Also the J\ — J2 model, in relation to Li 2 VOSi04 and Li 2 VOGe04 
compounds^, and the XY modeliS were analyzed with the same technique. All the studies performed have been based 
on a mean field decoupling scheme as the starting point to deal with the non-local interactions introduced by the JW 
transformation. In£ the mean field procedure was further supplemented by the inclusion of fluctuations in terms of 
an auxiliary gauge field with a leading Chern- Simons dynamics coupled to the lattice fermions. However, in spite of 
the partial success of the JW transformation, many problems remain open, in particular in connection to the study 
of frustrated systems such as the triangular lattice, etc. In some cases, the results obtained via a direct mean field 
treatment lead to results that are believed to be incorrect, like the appearance of a spin gap in the triangular lattice 
case (see the discussion m&) . The main problem associated with the JW approach is related to the implementation of 
the above mentioned mean field decoupling, which renders the description approximate. Another highly non-trivial 
problem is the construction of the lattice description of the Chern-Simons theory, which has been carefully studied 
for the square lattice case onljiii. 

It is the purpose of the present paper to propose a systematic self consistent mean field method for exploring the 
ground state (GS) of 2D lattice spin 1/2 systems, in a way that could be applied to arbitrary lattice topologies. The 
method can also be used in the presence of an external magnetic field, at finite temperature and even be applied to 
disordered systems. 
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II. JORDAN- WIGNER TRANSFORMATION IN 2-DIMENSIONS 



The Jordan- Wigner transformation in 2 spatial dimensions was originally proposed in^ as a generalization of the 
well known transformation in ID, and has been further developed inS^. It maps a set of spin 1/2 operators S p on 
lattice sites p into spinless fermion operators c p by 



c p exp 



»•£■ 



> c<c 



Sp 



cp exp 



CpCp 



1/2, 



(1) 



where = S™ ± iS' 2 ' are the usual spin raising and lowering operators and 9 qp is the argument of the vector drawn 
from site p to site q. The transformation is non-local, and sets a preferred quantization axes z. The spin operators 
l|T]l satisfy bosonic SU(2) commutation relations, while the Pauli principle ensures that they belong to the irreducible 
representation S = 1/2. Indeed, the only necessary ingredient that ensures the SU(2) commutation relations is the 
assignment of the phase factors which satisfies, for each pair of sites p, q 



= -1. 



(2) 



One should notice that there is a large freedom in choosing phase factors satisfying this condition (J2J). For instance, 
one could arbitrarily shift 9 pq — > 6 pq + 2kir with different integers k for each pair of lattice points p, q, or even perform 
an arbitrary simultaneous rotation for 9 pq and qp . Standard plane angles — tt < 8 < 7r measured from the x axis is 
just the simplest translation invariant choice on the flat infinite plane. It should be stressed that this large freedom 
does not alter the physical results, as long as all degrees of freedom are treated exactly. However, in any approximate 
treatment, this may introduce ambiguities that should be handled carefully, as we discuss below. 

One salient feature of the JW transformation is that no constraint is needed on the new variables (cf. for instance 
the Holstcin-Primakoff or Schwingcr bosons), but non- locality is the main stumbling block in the approach. 

The success of the JW transformation in 1 spatial dimension, in spite of being non-local, resides on the fact that XY 
nearest neighbors (NN) interactions become local in fermion variables; this is not the case in 2 dimensions. Indeed, 
consider the XY Hamiltonian on a given 2D lattice 



Hxy = J Y, ( S P S ^ + S P S t) ' 

<p,q> 



(3) 



where J is the exchange constant and the sum runs over all nearest neighbors < p,q > on the lattice. In terms of 
fermion variables the Hamiltonian reads 



Hxy = J (\ c l el&iP ' q)c i + H - c ) 



(4) 



<p,q> 



where 



a(p, q) = y^(g rg - Q rv )c\c r 



(5) 



(the ' on the sum indicates that 9 rr terms are absent). This phase is highly non local; in the ID case, this same 
expression becomes local due to the fact that the only two actual values for the angles are and tt. The non-locality 
in 2D is usually overtaken by the introduction of an auxiliary gauge field A^, which on the one hand represents the 
phases in eq.Q as the usual minimal coupling on the lattice, and on the other hand is governed by a Chcrn-Simons 
action. The Gauss law associated to the first order Chern-Simons action imposes a constraint which in anyon language 
attaches half a quantum flux to each fermion, providing the statistical transmutation of fcrmions into bosons. Then, 
a mean field treatment (known as average field approximation) of the gauge field can be done, leading in general to a 
quadratic NN interaction between fermions^. However, the Chern-Simons approach has serious difficulties when one 
deals with arbitrary lattice topologies (for example the triangular lattice) , and the associated mathematical problems 
are not yet solved. 
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We do not introduce such an auxiliary gauge field, but keep working with fcrmion variables. In order to perform 
numeric computations, one has to set a finite size lattice and impose suitable boundary conditions. We use periodic 
boundary conditions, this leading to a lattice on the torus. Moreover, the lattice size should be compatible with 
possible periodic configurations; in the case of a square lattice, size must be even in order to not interfere with the 
possible Neel order. 

Now, it is not straightforward to define the JW transformation on the torusi^, as the vector joining two different 
points is not unique . As one has to take care of condition @ , the vectors joining p with r and r with p must have 
arguments differing in tt. We have to choose a unique segment joining each pair of points p, r, and then draw both 
vectors along it. One can choose this segment by a criteria of minimal distance. However, there exist pairs of points 
on the torus that can be joined by two or more different segments with minimal distance and hence an ad-hoc criterion 
must be added. Any such criterion unavoidably breaks translation invariance, by preferring one segment over the 
rest. Naturally, we propose a criterion trying to minimize the violation of translation symmetry as follows: we set a 
principal finite size lattice and extend it on a plane by periodicity; for each point on the principal lattice we consider 
also its periodic copies. Now, given a pair of sites, we look for the shortest segment joining either the points or their 
copies; when such a segment is unique, the procedure is translationally invariant. For those pair of points where 
one can find more than one minimal distance segments, we choose the one with both ends belonging to the principal 
lattice, thus breaking translation invariance. Finally, the angles 9 pr and 6 rp are computed as the arguments of the 
vectors joining p and r along the chosen segment. For convenience we also define that 8 pp = 0, in order to handle the 
restriction on the sums in eqs.QISJ). 

As we mentioned in the Introduction, the JW transformation is exact but the resulting Hamiltonian is highly 
non-local and some kind of approximation is necessary to proceed. 

We propose here a variational approach to deal with the non-local phases in cq.Q) and the quartic terms that 
can arise from S z interactions. Working directly with fermion variables, we replace the local fermionic occupation 
numbers n p = c p c p by their expectation values in an arbitrarily chosen variational state. This procedure leads to a 
multi-parameter mean field approach, that will in turn be evaluated self-consistcntly. This is the subject of the next 
section. 



III. VARIATIONAL APPROACH, APPLIED TO THE XXZ MODEL 

To describe in full detail the method laid down above, we apply it to a generalized quantum spin 1/2 Heiscnbcrg 
antifcrromagnet in a square 2D periodic lattice, defined by the Hamiltonian 

Hxx z = JY1 ( S p S v + S P S ^ + AS p S D - h Y. S P' ( 6 ) 

<p,q> V 

where S p = (S p , S^, S p ) represents the S = 1/2 spin operator at site p, J > is the exchange constant and < A < oo 
the "XXZ" anisotropy parameter. The first sum in © runs over all nearest neighbors in the given lattice, while the 
last term represents the interaction with a transverse external magnetic field h. We work on a periodic rectangular 
lattice of size K = N x x N y . 

Using the JW transformation defined in eq.Q), the Hamiltonian can be written in terms of spinless fermions as 

h xxz = j \\i4 ei&Mc * + H - c -) + a (4 c p ~ \)] -hJ2(4 c p - \) > (7) 

<p,q> L P 

where the phase a(p,q) is defined in JSJ. Notice that the magnetic field h plays the role of a chemical potential for 
the JW fermions. In particular, we look for the ground state of the system (JJJ with fixed global magnetization M = 
(corresponding to h = 0). 

We implement a self consistent mean field solution by starting with a given fcrmion distribution profile {n p }t (which 
can be random or guided by some ansatz) on the lattice, 

(n p ) = rip, (8) 

which has to satisfy a global constraint to provide the given magnetization (here J2 n p = K/2 corresponds to M = 0). 
We then replace the operator a(p) by its expectation value 



?)) = ^2(°rq - rp ) TV, 
r 



(9) 



4 



where the angles 9 pq are assigned following the criterion presented in the previous section. To be precise, the principal 
lattice can be defined by indexing each site by a position pair and setting the range i = - ■ ■ N x — 1, j = 

■ • • N y — 1. Periodic boundary conditions are then expressed by = (i + N x ,j) = (i,j + N y ). 
Regarding the Ising term 

SpSq = C p C p C q C q ~ 7^\) C V — ~^q C 1 ~ 4 (^) 

in eq.JJJl, it is quartic in fermion operators, so we also treat it in mean field. In order to approximate the first term 
in I|1U|) with a quadratic expression we propose the following 

supported by best results in a posteriori evaluation of the GS energy (some other possibilities arc discussed in-). 
At this step, the Hamiltonian can be written as 

JMF) 



H { xxJ({n p }) = J2 clJ pq ({n p })c q + C (12) 



where 



L e "*(p,<i) if < p ; g > nearest neighbors 

Jp* = { ^(E„ eifffc6ors q (n q ) - 4) if p - q (13) 

otherwise 

and C = KJA/2. 

The main idea of the present paper is to provide a systematic way to compute an approximation to the true GS. 
We first find the GS for the quadratic H^fJ ({%>}) by solving the one particle (IP) spectrum and filling the system 
with the lowest energy IP states, up to the proper filling fixed by the total magnetization M. Then we compute from 
this approximate GS a new set of local densities n' p = (GS\c p c p \GS), which we use as a new input in l|12|) and iterate 
this procedure looking for a fixed point configuration for the density profile, i.e. a set of local densities {n*} satisfying 



n' p ({n* q }) = n* p . (14) 

The existence of a fixed point solution for this mapping and its eventual dependence on a given initial configuration 
is not at all obvious and has to be studied numerically. 

In order to proceed with the method, H^fJ^rip}) can be written in diagonal form 

K 

Hxxz({ n p}) = e (^) d\.dk + constant (15) 
fc=i 

where arc the IP eigenvalues of the quadratic part of H^j^J . Notice that k is just an integer index over the 
spectrum, not to be confused with the lattice momentum. Moreover, we order the eigenvalues asccndcntly. 
The operators du are related to c p by 

dk = Y^Qlp c P ( 16 ) 

p 

where Q p u is the matrix of eigenvectors of J vq . We compute both and Q p k numerically. Being Q unitary, the set 
of dk operators satisfy fermion commutation relations, {dfc,dj,/} = Skk'- Moreover, the total fermion number operator 
satisfies 

K K 
p=l k=l 

making it easy to control the filling in terms of the new fermions. 
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We now construct the approximation to the quantum GS as the half-filled state that minimizes the energy, namely 

K/2 

\GS) = l[dl\0) (18) 
fe=i 

Notice that this is a well defined quantum state of K/2 particles, except for casual degeneracy of the IP spectrum at 
the Fermi level . This is not the case for the XX Z model on the square lattice (see details below). 
From \GS) it is now easy to compute the approximate GS energy as 

E GS = (GS\H^\GS) = J2 ^ + C - ( 19 ) 

k<K/2 

Also the local occupation numbers can be computed in this approximate GS as 

n' p = (GS\clc p \GS)= J2 %kQ P k- (20) 

k<K/2 

With these occupation numbers we start again the procedure: compute J pq in MF, diagonalize the new H^^fJ , etc. 

We have found after thorough numerical investigations that a fixed point solution for ea. (|14|l always exists, but 
metastable solutions can also appear, depending on the initial configuration one chooses. In any case, one can 
distinguish metastable solutions from the best GS approximation simply by comparing their energies. Moreover, we 
describe below how this drawback can be naturally solved by introducing a thermal bath to kick the system out from 
the vicinity of metastable states. 

Indeed, one can consider the effects of finite temperature by replacing the proposed ground state l|18fl by a thermal 
state l^), compatible with the Fermi-Dirac IP energy distribution at a given temperature, 



"< e > " l + c*p W . -.-)) ' <21) 

where e is the IP Fermi energy at half filling, and (3 = 1/kT is the inverse temperature. In detail, this thermal state 
|<3/ p) is constructed as 



i*/3>=n4i°>> ( 22 ) 

fcGcr 

where a is a set of K/2 IP states chosen with probability n(e(k)) from some random simulation. 

An exploration of the Hilbert space of the system by constructing a thermal state from a starting fermion dis- 
tribution, computing from it the new local fermion distribution and again constructing a thermal state should be 
considered as a thermalization at the given temperature. It provides a source of thermal noise that has proven to 
help the system in finding lower energy fixed points. 

The thermalization can be done through several steps at a given temperature, and then quenching to the pure 
quantum regime (T = 0), or it can be implemented by gradually lowering T (annealing). 

Besides, results at finite T can also be achieved, by constructing an statistical ensemble of microscopic states 
compatible with T. Observables should then be computed as averages over the statistical ensemble. We do not 
attempt to complete this program in the present paper. 



IV. RESULTS 



We have tested the iterative approach described in the previous section with the well known anisotropic XX Z 
model on periodic 2D square lattices of size up to 20 x 20 sites, at zero total magnetization. The sizes of the 
lattice that we explored are by no means an upper limit, as our computations were made on a modest computer. 
The anisotropy parameter A has been explored in a range from 0.05 to 1.5, including the isotropic SU(2) case 
(A = 1, Heisenberg model). As starting configurations {n p } we have used random, uniform, and different amplitude 
staggered distributions. We performed several iterations and analyzed the evolution of the local fermion profile and 
the approximate GS energy. We report the results in terms of spin variables, noting that the local fermion occupation 
represents the local magnetization as m z (p) = (n p — \)- 
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Working at T = 0, we have found that in general, from different starting configurations, the system rapidly finds a 
Neel order as stable ground state approximation, after 15 ~ 20 iterations. The Neel order parameter, usually defined 
as the staggered or sublattice magnetization m z , depends on the anisotropy parameter A. Fluctuations around this 
staggered magnetization are typically of order 10~ 8 . In figure Q we plot the Neel order parameter m z of the fixed 
point solution for different values of A, for several lattice sizes. Finite size effects are noticeable for lower values of A, 
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FIG. 1: 

shown. 



Neel order parameter at fixed points, in function of the anisotropy A. Several lattice sizes and the scaling limit are 



so we also show the results of a finite size scaling m 2 (oo) of our data, fitted with a power law m z (K) = m z (oo) + c/K a . 
The corresponding GS energies per site are shown in figure (J2J) where one observes that scaling with the system size 
is clearly less important. We have observed that the IP spectrum of the mean field Hamiltonian (|12fl presents a gap 
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FIG. 2: GS energy, in function of the anisotropy A. Data corresponds to configurations plotted in figure 

2m z J A for Neel ordered configurations, at the half filling Fermi level. This is in agreement with (— ) and makes the 
construction of \GS) in ca. l|18(l unambiguous. 

In the case of random initial distributions, metastable configurations can show up; a detailed inspection of the 
local magnetization in these cases reveals the formation of antiferromagnetic domains, that is the presence of the two 
possible Neel configurations in different regions. In figure Q we show an example of such domains, at two different 
stages of a sample evolution. It is natural to expect that larger lattices favor the formation of these domains, as 
it indeed is observed. These configurations have higher energy than the uniform Neel state and correspond then to 
metastable configurations; correspondingly, they are not presented in figures QI^Jl. 
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FIG. 3: Occupation patterns for a metastable configuration, where antiferromagnetic domains appear. The size of the points 
is proportional to the local fermion occupation number. This configurations occurred on a 20 x 20 lattice, with A = 1.1, after 
10 steps of thermalization at T = 0.2 J, and 10 (left panel) or 20 (right panel) more steps of GS search at T = 0. The smaller 
domain is seen to decrease in size under the simulated evolution. 



When a thermal bath is simulated on random initial configurations, we have observed that metastable configurations 
are less likely to appear. After thermalization we let the system to cool down by either quenching or annealing as 
described in Section III, and complete the iterations at T = 0. In fact, a few steps (~ 10) of thermalization with 
sufficiently high T completely avoid domain formation and lead to a unique fixed point mean field configuration; the 
required temperature is higher for larger lattices, being of the order of J for the lattice of 20 x 20. We have checked 
that under general circumstances, quenching provides the fastest convergence method to the minimum energy state. 
An example of the evolution of the Neel order parameter from an initial random configuration, under thermalization 
with different temperatures, is shown in figure 
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FIG. 4: Example of the evolution of the Neel order parameter from a sample initial random configuration. Data corresponds to 
the system depicted in fig.|3"). with a vertical line separating the thermal evolution and the T = evolution. Error bars indicate 
the standard deviation of local magnetization from Neel order (reduced by a factor of 5 for clarity). Insufficient thermalization 
can lead to metastable configurations or to very slow convergence, while higher temperature dramatically improves convergence 
towards an ordered configuration. 

The results of the present MF computation show all the features expected for the Hcisenberg antiferromagnet on 
the square lattice. They are of course not comparable to accurate numerical techniques^, but are in qualitative 
agreement with results from previous studies. In particular, in the scaling limit we obtain no Neel order for small 
anisotropy A, where the system presumably has XY order. We can estimate a critical value A* « 0.2, above which 
Neel order develops. For the isotropic Heisenberg point A = 1 we obtain a sublattice magnetization m z = 0.3453, 
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with ground state energy per site Eqs/K = —0.5683 J, to be compared for instance with corresponding Quantum 
Monte Carlo values of 0.307 and -0.6694Jil 



V. CONCLUSIONS 



We have presented a self consistent MF procedure for exploring the quantum ground state of any S = 1/2 spin system 
on a 2D lattice. When tested on the XXZ model on a square lattice, the method provides the correct qualitative 
description of the system, with no a priori ansatz for any kind of order. We computed the values for the sublattice 
magnetization and GS energy for a wide range of values of the anisotropy parameter, which compare qualitatively 
well with the available numerical data, at least for A = 1 where most accurate data is available. Moreover, we have 
found that the sublattice magnetization as a function of the XXZ anisotropy shows the correct qualitative behaviour, 
expected from a spin wave analysis^. 

The present approach has a more general scope than previous MF computations, in the sense that it can be applied 
to any lattice topology, irrespectively of the appearance of frustrating units, a fact that prevents the applicability of 
one of the most powerful numerical techniques such as Quantum Monte Carlo. A magnetic field can be trivially added 
as a chemical potential for the JW fermions and hence magnetization curves could be obtained. Since the method 
is not based on any periodicity of couplings, it can be well suited to study disordered quantum spin systems, at the 
only price of increasing the CPU time. Last but not least, the approach is naturally well suited for the study of the 
thermodynamics of these systems, since temperature can be added in a simple way. 

Among other situations, it would be interesting to apply this technique to the Heisenberg quantum AF on the 
triangular lattice, where there is disagreement between Chern-Simons MF predictions^ and numerical data about a 
magnetization plateaux at zero magnetization. Another case of interest the kagome lattice, where a quantum spin 
liquid is believed to be realized^ (see alsoi&). This issues will be investigated elsewhere. 
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